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We present a novel technique well suited to study the ground state of inhomogeneous fermionic 
matter in a wide range of different systems. The system is described using a Fermionic Shadow 
wavefunction (FSWF) and the energy is computed by means of the Variational Monte Carlo tech- 
nique. The general form of FSWF is useful to describe many-body systems with the coexistence 
of different phases as well in the presence of defects or impurities, but it requires overcoming a 
significant sign problem. As an application, we studied the energy to activate vacancies in solid 
^He. 



The microscopic theoretical description of inhomoge- 
neous Fermionic systems is a long-standing chaUenge. 
Among such systems we can include defective quantum 
crystals (e.g. ^He crystals or electron Wigner crystals 
in presence of vacancies or defects) or where an ordered 
and a disordered phase coexist like, for instance, fluid 
and crystal. The main difficulty consists in the fact that 
one has to deal with a wavefunction that combines the 
antisymmetry required by the Pauli principle and the in- 
homogeneity of the system itself. 

While mean field methods are very efficient in dealing 
with homogeneous phases (e.g. an extensive and perfect 
solid), the phase coexistence or the description of local 
defects presents difficulties. The main reason is the local- 
ity of the inhomogeneity. Therefore, an explicit descrip- 
tion of the wavefunction seems a much better approach. 

The rigorous microscopic evaluation of the vacancy for- 
mation energy in ^He is one of the problems that most 
suffers from the limitations of standard theoretical tools. 

The problem was successfully solved for many-Boson 
systems several years ago, by means of the so-called 
Shadow Wave Functions (SWF)[l|,[li], a class of wave- 
functions based on the introduction of auxiliary degrees 
of freedom, which was successfully applied to a variety of 
inhomogeneous phases of "^Hejl, 0] and P-H2 0: 
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where R = {ri, tn} are the coordinates of the TV con- 
stituents of the system, and S = {si, sn} are auxiliary 
degrees of freedom called "shadows" . (f)p and (j)s are two- 
body correlation factors (of the so-called Jastrow form) 
for particles and shadows, respectively, and S(R, S) is 
the kernel describing the correlations between particles 
and shadows. Here, as in most applications, we take the 
kernel to be a Gaussian: 

S(R,S) = exp(-c(R- S)2). (2) 
The main properties of the SWF are i) the fact that 



it introduces correlations to all orders via the integra- 
tion over the auxiliary degrees of freedom and ii) the 
fact that despite its manifest translational invariance, it 
can describe phases in which the translational symme- 
try is broken (solids, interfaces, defects). However, the 
extension to many-Fermion systems is hard. A straight- 
forward extension (which we term ASWF) based on the 
antisymmetrization of the particle degrees of freedom was 
proposed several years ago for the study of homogeneous 
■^HeJ,], and later apphed to the homogeneous electron 
gas[ 



i^ASWFi^) = n A[0fc(rz)]0p(R) 

X y"exp(-c(R- S)2)0,(S)dS, (3) 



where Di[<j)k{vi)\ is a Slater determinant of particle or- 
bitals. 

This form has an evident drawback. Once the orbitals 
are specified, as e.g. plane waves satisfying the Born-von 
Karman conditions for an extensive system, the nodal 
structure remains unchanged even if the auxiliary degrees 
of freedom provide the correlations necessary to break the 
symmetry so as to correctly describe a set of localized 
particles. 

The only viable solution to the problem is to introduce 
an antisymmetric wavefunction in which crystallization 
can be described without explicit symmetry breaking, 
but that at the same time can develop a correct nodal 
structure according to the phase described. This goal 
can be achieved by writing a Shadow Wave Function in 
which the antisymmetry is imposed on the auxiliary de- 
grees of freedom, therefore maintaining explicit correla- 
tions among the particles only in the symmetric part of 
the function. The so-called Fermion-Shadow Wave Func- 
tion (FSWF) assumes the following form: 
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fswf(R) = </'p(R) J exp(-c(R-S)2) 

X H Di[Ms,)](f>s{S)dS. (4) 



It is possible to prove that ipFSWF is antisymmetric un- 
der the exchange of two particles of like spin 0. The 
main difference between ASWF and FSWF comes from 
the fact that the latter develops a nodal structure for 
the particle degrees of freedom that depends on the in- 
tegration over the shadow degrees of freedom, and in- 
cludes effects of correlations to all orders. In particular 
it is easily proved that when particles and shadows are 
strongly localized by the effect of the two-body correla- 
tions among the shadows themselves, the wavefunction is 
closely approximated by a determinant of Gaussians con- 
necting each particle to each shadow, which is obviously 
closer to the structure of the wavefunction expected for 
a quantum crystal. 

There is a very high technical price to pay in or- 
der to exploit FSWF in computations. In a Variational 
Monte Carlo (VMC) calculation, the absolute square of 
the wavefunction is used as a probability density, P(R), 
and the local energy is averaged over the sampled config- 
urations. Thus 

JdIir{^Mn)EL ^ J dIlP{IL)EL 
/dRV'*(R)V(R) /dRP(R) ' ^ ' 

where El — ■ip~^Hip is the local energy of the system. 
The integral is evaluated by generating configurations ac- 
cording to P — that are sampled using the Metropo- 
lis algorithm. 

In using SWF, one constructs ■(/'*(R)'0(R) by integrat- 
ing over two sets of shadow variables, S and S'. Define 
the integrand of SWF as 

Z(R,S,S') = 0^(R) 

X exp(-c(R - Sf - c(R - S')^)(/),(S)0,(S') (6) 

and choose some probability density function, 
P(R, S,S'). Sampling in the usual way, a generic 
operator can be computed as 

_ / dndSdS'PCR, S, S')wCR, S, S')O(R) 



where 
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A reasonable choice for ordinary SWF is the integrand 
itself: 



with w = 1 and P = P = Z . Similarly in 4'aswf the 
integrand is positive, and can be used for P(R, S, S'). 

When using V'fsivf, the integrand is not positive defi- 
nite, and sampling it is not possible. It is always possible, 
however, to sample a suitable probability distribution 
and compute a weighted average. The most straight- 
forward choice in this case is the absolute value of the 
integrand. Let 

Q(R, S, S') = (/.2(R) exp(-c(R - S)^ - c(R - S')') 

X n DiiM^,)] n Di[Ms'^]MS)US'), (10) 



1=1,1 i' = T,i 

P(R,S,S') = |Q(R,S,S')|, 
-(^'"'")-|Q(R,S,S')r^'- 



(11) 
(12) 



It should be noted that the normalization integral, while 
containing positive and negative terms, is always posi- 
tive by construction. However, the speed of convergence 
of the integral strongly depends on the fluctuations in 
sign of the integrand. In particular, the intrinsic vari- 
ance of the integrand might become excessively large, and 
prevent the computation of an average with acceptable 
statistical errors. This is particularly true for disordered 
systems, such as liquids or disordered solids, where the 
wavefunction has strong variations in space. 

A simple reorganization of the calculation produces a 
dramatic improvement in the Monte Carlo efhciency. In 
applying the Metropolis method to shadow wavefunction, 
including the ASWF variant, the usual procedure is to 
sample new values of R, S, and S' in turn. It is always 
true that the integrals in Eq. ([7]) over S and S' for fixed 
R are positive. This suggests that a change in the or- 
der of summation might be useful for the FSWF class of 
functions, especially with disorder, by propagating the 
shadows S and S' for M steps (with M big enough) for 
fixed R. That is, we expect that the sum of the ±1 weight 
of M steps. 



WsWs' 
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P(R,S,S') ^ Z(R,S,S'), 



(9) 



will be usually positive. In fact, increasing M in the 
more difficult cases where exchanges of sign often appear, 
gives weights WsVFs' usually positive and significantly 
different from zero. The algorithm becomes: i) sample a 
configuration R of particles, ii) sample M configurations 
of S, iii) sample M configurations of S', iv) combine all 
the weight factors and accumulate the local energy for 
the average and variance, iterate from i) to iv) until the 
convergence is reached and the variance is low as desired. 
For a calculation of the crystalline phase with no empty 
sites we typically sampled configurations of particles, and 
using M = 1 or M = 100 does not significantly change 
the result. If a vacancy is present in the system and 
the sign frequently changes, we usually sampled 5 x 10^ 
configurations and for each one M = 1.5 x 10^. 
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As an illustration of the capabilities of FSWF, we stud- 
ied the ground state of solid ^He in the range of molar 
volumes between 20 and 24 cc/mol. A comparison with 
previous variational estimates based on standard anti- 
symmetric wavefunctions is given below. As a next step, 
we studied the vacancy formation energy by computing 
the energy in the presence of an empty site. As already 
pointed out, standard imaginary time projection calcu- 
lations require that a wavefunction imposing crystalliza- 
tion be used, preventing one from taking properly into 
account all the phenomenology related to lattice relax- 
ation and vacancy mobility. As model He-He interaction 
we chose the Hartree-Fock dispersion HFDHE2 potential 
by Aziz et al. [lOl] . which gives an overall description 
of the equation of state in good agreement with experi- 
ments, though it does not introduce explicit three-body 
terms. In the Jastrow functions ips the pseudopotential 
Usifij) was taken as the rescaled particle-pair potential, 
v{rij) namely Us{rij) = av[(3rij) {a and (3 are additional 
variational parameters), while in the (j)p we used a McMil- 
lan form lZII combined with a summation over a basis as 
in Ref. All the variational parameters entering in 

the wavefunction were optimized at each density using 
an energy-variance minimization technique due to C.J. 
Umrigar and M.P Nightingale applied to the system with 
no defects. 





E/N(54) 


T/N(54) 


E/N(53) T/N(53) 


0.419 


0.422(1) 


23.947(1) 


0.69(1) 


23.71(3) 


0.427 


0.548(2) 24.608(1) 


0.78(8) 


24.6(2) 


0.438 


0.955(1) 


26.005(2) 


1.26(1) 


25.76(3) 


0.457 


1.556(1) 


27.986(2) 


1.844(8) 


28.05(2) 


0.479 


2.455(1) 


30.482(2) 


2.801(7) 


30.60(2) 


0.503 


3.481(1) 


32.487(2) 


4.127(7) 


32.35(2) 



TABLE I: Total and kinetic energy per particle at different 
densities in the solid b.c.c. phase for the crystal with no 
defects (54 atoms) and with the presence of an empty site (53 
atoms). All the energies are expressed in K. 

We report the energy of 54 atoms in table HI The en- 
ergy per particle is also displayed in Fig. [T] where we 
compared our results (circles) with those found in Ref. 
computed using ASWF starting from a normal anti- 
ferromagnetic order NAF (diamonds), and including ex- 
changes (triangles), and with the more accurate Diffusion 
Monte Carlo (DMC) results of Ref. [ll| (squares). As 
can be seen FSWF provides the lowest of the variational 
estimates of the energy. The DMC energies are lower at 
each density by a constant value of about IK. 

The vacancy formation energy at constant pressure for 
a system with N particles at a fixed density p can be 
computed as 




0,42 0,43 



Ai;„ = E{N-l,Ni = iV)- 



N -I 
N 



EiN, Ni=N), (14) 



FIG. 1: (color online) The FSWF energy per particle (black 
circles) as a function of the density. The result is compared 
with the two results provided by ASWF of Ref. 7] (blue 
trian gles and red diamonds) and with the DMC results of 
Ref. 111. 



where the number of lattice sites Ni is conserved and 
the density of the two systems is the same. The vacancy 
formation energy includes contributions from lattice re- 
laxation and tunneling that cannot be accounted for by a 
wavefunction with an underlying lattice structure (such 
as a Jastrow-Nosanow wavefunction). The computation 
of the energy for the system with — 1 particles is per- 
formed by removing one particle and one shadow from 
the trial wavefunction. It has to be noted that dropping 
one shadow from the shadow determinant means hav- 
ing a hole state in one of the determinants of Eq. H) 
In the case of an open-shell configuration where one or 
more single particle states are not filled one needs to per- 
form the calculation using twist average boundary con- 
ditions [3) EBl , so that the total wavefunction has zero 
total momentum. Using FSWF, we conjecture that it is 
not important which particular shadow orbital is omit- 
ted, because the total momentum of the system is always 
conserved. 

We tested this by repeating the same calculation by 
removing different single-shadow states from the kernel. 
In particular the energy of the system where orbitals 
with different vector number n were removed is the same 
within statistical error. We stress that, by contrast, in 
using a normal many-body wavefunction with no shad- 
ows or using ASWF the energy would depend upon the 
unfilled single-particle orbital. We assume for now that 
the system containing one empty site is well described by 
the same wavefunction of the system with the complete 
crystal. The modified structure of the system with a va- 
cancy is modeled by the shadow extra variables. There- 
fore for 53 atoms we used the same parametrization of 
(j)p, (j)s and the coupling constant c entering in Gaussians 
of the system with 54 atoms. 

The energies of the system with a vacancy are reported 
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FIG. 2: (color online) The vacancy formation energy Ai?« as 
a function of the density using FSWF (black points) and by 
keeping the shadows fixed. Some experimental data from Ref. 
is also reported for comparison. See the text for details. 



in table H] The vacancy formation energy obtained using 
Eq. [M] is given in Fig. [2l where we included the same 
calculation using the static shadows (red points). In the 
latter case the shadow degrees of freedom are kept fixed 
on the lattice sites so their effect is switched off. This 
corresponds to using a variational wavefunction of the 
antisymmetric Jastrow-Nosanow type. FSWF is more ef- 
fective as is clear from the figure. The vacancy formation 
energy computed by means of FSWF is larger than the 
experimental data (blue points) taken from Ref. [13] (see 
also Ref. pj[ and references therein). The discrepancy 
can be attributed to several possible sources. First of all 
the calculation might be affected by strong finite-size ef- 
fects. In fact the effective concentration of vacancies in 
the system is rather high (i.e. 1/N), and this might imply 
a contribution to the vacancy-formation energy coming 
from a vacancy-vacancy interaction. There is also ad- 
ditional room for improvement in the overall variational 
description. For example, the parameters could be re- 
optimized in the presence of the vacancy, or a more so- 
phisticated version of the wavefunction including a local- 
density dependence of the two-body correlations might 
be usedQ. 

In conclusion, in this Letter we present a novel varia- 
tional wavefunction to study fermionic systems with im- 
purities. We describe the Fermionic Shadow wavefunc- 
tion that we used to compute the equation of state of solid 
^He in the b.c.c. phase, and the vacancy formation en- 
ergy as a function of the density. We stress the fact that 
using standard wavefunctions it is not possible to cor- 
rectly study systems with the presence of defects or im- 
purities like a vacancy. In particular the theoretical study 
of ^He with vacancies requires a correct description of re- 
laxation and tunneling effects that cannot be addressed 
by using standard forms of wavefunctions. Within the 
variational framework we also computed the pair distri- 



bution functions between atoms and the density around 
a vacancy, and the calculation of other properties is pos- 
sible, but this, as well as the technical improvement of 
our present method, will be the subject of future work. 
Using FSWF it is possible to study '^He with the pres- 
ence of impurities of ^He as well as the mixture of the two 
gases, and to move near the region where the solid and 
liquid phases start to coexist. Work in these directions 
is in progress. 
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